### Needed packages
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggmap)
library(scatterpie)

### Needed variables

# Make a nice color pallete and legend order for all plots

##  health status
my_cols=c("darkgreen", 
              "darkred", 
              "orangered1", 
              "cadetblue", 
              "tan", 
              "beige", 
            #  "burlywood4", 
              "coral", 
              "aquamarine3", 
              "gray70", 
              "black")

desired_order=c("healthy", 
                "ozone", 
                "ozone_and_other", 
                "others_combined", 
                "drougth", 
                "fungi", 
             #   "insect", 
                "worm", 
                "acid_rain", 
                "other", 
                "dead")
  
 spanish_labels=c("Sano", 
                  "Ozono", 
                  "Ozono y otros", 
                  "Otros combinados no-ozono", 
                  "Sequía", 
                  "Hongos", 
               #   "Insectos", 
                  "Gusano de seda", 
                  "Lluvia acida", 
                  "Otro", 
                  "Muerto")
 
## Damage percentage
my_cols2<-c("gold2", "chocolate1", "orangered", "red4", "darkorchid4")
 
desired_order_percentage<-c("less than 10%", "10 to 40%", "40 to 50%", "50 to 70%", "more than 70%")

Este reporte explora los resultados del monitoreo participativo del estado de salud de árboles de oyamel en el Parque Nacional Desierto de los Leones y sus zonas de influencia en Santa Rosa Xochiac.

Los datos corresponden al resultado del muestreo participativo realizado por 12 brigadistas de Santa Rosa Xochiac, utilizando kobo-conabio como parte del proyecto 308488 Monitoreo y manejo para la conservación de bosques aledaños a la CDMX afectados por contaminación atmosférica de la convocatoria FORDECYT 2019-5.

Los datos del muestreo corresponden a los datos colectados con kobo y limpiados previamente con el script 1_preprocesamiento_datos_kobo.Rmd.

# load data
muestreo_tidy<-read.delim("../data/kobo/muestreo_dic2020_tidy.txt", header = TRUE)
parcelas_tidy<-read.delim("../data/kobo/parcelas_dic2020_tidy.txt", header = TRUE)

# pivot long parcelas data to have health data as a single variable
parcelas_long<-pivot_longer(parcelas_tidy, 
                            cols = healthy:worm, 
                            names_to = "tree_health_simplified",
                            values_to = "n_trees")

Los datos analizados aquí corresponden sólo a los árboles que fueron aprovados durante la validación revisando manualmente las fotografías en kobotoolbox. Del total de 1771 árboles muestreados, 1718 fueron aprovados en la validación.

muestreo_tidy<- filter(muestreo_tidy, X_validation_status=="validation_status_approved")

Distribución del estado de salud de los árboles por parcela

La siguiente figura muestra el total de árboles muestreados en cada parcela de 10x10 m, y cuántos de estos están bajo alguna categoría de daño:

p <- ggplot(parcelas_long, aes(x=plot, y=n_trees,     fill=tree_health_simplified)) +
  geom_bar(stat="identity") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") 
  

p + theme_bw() +
  ggtitle("Estado de salud de los árboles por parcela") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  labs(x="Parcelas", y= "número de árboles") +
  theme(text = element_text(size = 25))

Esta es la distribución de las 48 parcelas:

# code adapted from https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-google-map.html

## configure google api

# You first need to register your api key in https://cloud.google.com/maps-platform/#get-started and follow instructions. The geocoding API is a free service, but you nevertheless need to associate a credit card with the account. Please note that the Google Maps API is not a free service. There is a free allowance of 40,000 calls to the geocoding API per month, and beyond that calls are $0.005 each.
# after you obtain your api, save it in /scripts/api_key.api (not shown in this repo por obvious reasons).

# if you get the following error when running get_map():

#"Error in aperm.default(map, c(2, 1, 3)) : 
#  invalid first argument, must be an array " 

# check this troubleshooting: https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-google-map.html

##  load and register api
api <- readLines("api_key.api")
register_google(key = api)

## plot map
# get map
sat_map = get_map(location = c(lon = -99.3060, lat = 19.2909), zoom = 14, maptype = 'satellite', source = "google")

# plot sampled plots
p_satmap <-  ggmap(sat_map)
p_satmap + geom_point(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude),
                      color="red") +
          geom_text(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude,
                          label=plot),
                      color="white",
                 #    check_overlap = TRUE,
                      hjust = 0, vjust=1, nudge_x = 0.0005,
                 size= 5) +
    theme(text = element_text(size = 30))

Según nuestro muestreo, el daño por ozono se distribuye espacialmente de esta forma:

# plot pies in map
p_satmap +
geom_scatterpie(data=parcelas_tidy,
                aes(x=X_coordinates_longitude,
                    y=X_coordinates_latitude,
                    group=plot),
                pie_scale = 1.5,
                cols=desired_order,
                color=NA,
                alpha=1)  +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") +
  theme(text = element_text(size = 20))

Haciendo un acercamiento a la zona central:

# zoom plot
sat_map = get_map(location = c(lon = -99.303, lat = 19.2890), zoom = 16, maptype = 'satellite', source = "google")

# plot
p_satmap <-  ggmap(sat_map)

# plot pies in map
p_satmap +
geom_scatterpie(data=parcelas_tidy,
                aes(x=X_coordinates_longitude,
                    y=X_coordinates_latitude,
                    group=plot),
                pie_scale = .7,
                cols=desired_order,
                color=NA,
                alpha=1)  +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") +
  theme(text = element_text(size = 30)) +
  
  # add plots id
geom_text(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude,
                          label=plot),
                      color="white",
                     check_overlap = TRUE,
                      hjust = 0, vjust=1, nudge_x = 0.0005,
                 size= 5)

En total se muestrearon 1718 árboles de los cuales 1676 se encontraban vivos. De estos, el 7.4% presentó daño por ozono, y el 29.36% presentó daño por ozono en combinación con otro tipo de daño, lo que hace al ozono la fuente de daño más abundante de estos bosques.

A continuación analizamos el daño por ozono en términos de porcentaje de árboles dañados por parcela.

Según la altitud:

# Create new variable with porcentage of ozonoe damage
parcelas_tidy<-parcelas_tidy %>% rowwise() %>% 
                     mutate(., 
                      total=sum(healthy,ozone,ozone_and_other,
                          drougth, acid_rain, other,
                          others_combined, dead, fungi,
                          # insect, 
                          worm)) %>%
                    mutate(perc.ozone= sum(ozone, ozone_and_other)/total)

#plot
p <- ggplot(parcelas_tidy) +
     geom_point(aes(x=X_coordinates_altitude,
             y=perc.ozone))

p + theme_bw() +
  ggtitle("Daño por ozono según altitud") + 
  labs(x="Altitud de la parcela", 
       y= "Porcentaje de árboles con daño por ozono")+
       theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  theme(text = element_text(size = 16))

O según la latitud. Esto es relevante porque latitudes más al sur están más lejos de la CDMX y por ende de la fuente de contaminantes:

p <- ggplot(parcelas_tidy) +
     geom_point(aes(x=X_coordinates_latitude,
             y=perc.ozone))

p + theme_bw() +
  ggtitle("Daño por ozono según latitud") + 
  labs(x="Latitud de la parcela", 
       y= "Porcentaje de árboles con daño por ozono")+
       theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  theme(text = element_text(size = 16))

Sin embargo, ni la altitud ni la latitud parecen influir en la distribución del daño por ozono.

Distribución del estado de salud de árboles individuales

Examinemos el daño por ozono dependiendo de si la planta fue reforestada o no, y de si se encuentra cubierta o expuesta.

p <- ggplot(muestreo_tidy, aes(x=tree_exposition, 
                      fill=tree_health_simplified)) +
       geom_bar(stat="count") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +  
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
  scale_x_discrete(breaks=c("cover", "exposed"),
        labels=c("cubierto", "expuesto"))
p + theme_bw() +
  ggtitle("Estado de salud de los árboles") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  labs(x="", y= "Número de árboles") +
  theme(text = element_text(size = 20))

En términos porcentuales:

p <- ggplot(muestreo_tidy, aes(x=tree_exposition, 
                      fill=tree_health_simplified)) +
       geom_bar(stat="count", position = "fill") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +  
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
  scale_x_discrete(breaks=c("cover", "exposed"),
        labels=c("cubierto", "expuesto"))
p + theme_bw() +
  ggtitle("Estado de salud de los árboles") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  labs(x="", y= "Porcentaje de árboles") +
  theme(text = element_text(size = 20))

Examinemos el daño según la altura de los árboles:

p <- ggplot(muestreo_tidy) +
 scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
theme_bw()

p2<-p + geom_histogram(aes(x=tree_heigth, 
                      fill=tree_health_simplified))  +
    labs(x="Altura del árbol (m)", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de los árboles según su altura") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) 
p2

p2 + 
  facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))

Ahora nos enfocamos en los árboles <15 m porque esta es la altura máxima en la cual podemos contar los nodos (lo que sirve para estimar la edad, como se explica más adelante):

p <- filter(muestreo_tidy, tree_heigth<15, tree_nodes>0) %>% 
     ggplot(.) +
     scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
theme_bw()
p + geom_histogram(aes(x=tree_heigth, 
                      fill=tree_health_simplified))  +
    labs(x="Altura del árbol (m)", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de árboles < 15 m según su altura") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold"))

Cada nodo se forma en un año de crecimiento, lo que permite estimar la edad de los árboles. Por lo tanto, la figura de abajo nos dice cómo se distribuye el daño en diferentes edades de árboles:

pnodos <- p + geom_histogram(aes(x=tree_nodes, 
                      fill=tree_health_simplified))  +
    labs(x="Número de nodos", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de árboles < 15 m \n según su número de nodos") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold"))
pnodos

Exploremos la figura anterior pero dividiendo en si fueron reforestadas o no:

pnodos + facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))

Al examinar esto en términos porcentuales, podemos observar que la cantidad de árboles con daño por ozono (o daño por ozono más otros daños) aumenta con la edad:

pnodos <- p + geom_histogram(aes(x=tree_nodes, 
                      fill=tree_health_simplified),
                      position= "fill", binwidth=1)  +
    labs(x="Número de nodos", y= "Porcentaje de árboles") +
    theme(text = element_text(size = 20)) +
    ggtitle("Estado de salud de árboles < 15 m \n según su número de nodos") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold"))
pnodos

pnodos + facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))
## Warning: Removed 99 rows containing missing values (geom_bar).

También podemos explorar la distribución del daño por categorías diamétricas, donde encontramos un patrón parecido al anterior.

Por número de árboles:

# changer order of levels of tree_diameter_category
muestreo_tidy$tree_diameter_category<- factor(muestreo_tidy$tree_diameter_category,
levels=c("", "0.5_cm", "1_cm",   "2_cm", "6_cm", "10_cm", "30_cm" , "40_cm")) 

levels(muestreo_tidy$tree_diameter_category)<-c(NA, "0.5_cm", "1_cm",   "2_cm", "6_cm", "10_cm", "30_cm" , "40_cm")

# subset data and plot
filter(muestreo_tidy, tree_heigth<15, !is.na(tree_diameter_category)) %>% 
     ggplot(.) +
     geom_bar(aes(x=tree_diameter_category, 
                      fill=tree_health_simplified),
                      stat="count")  +
    labs(x="Categoría diamétrica", y= "Número de árboles") +
    theme(text = element_text(size = 15)) +
    ggtitle("Estado de salud de los árboles según su diámetro") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) + 
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
    facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
    theme(axis.text.x = element_text(angle=90)) + theme_bw()

En términos porcentules:

# subset data and plot
filter(muestreo_tidy, tree_heigth<15, !is.na(tree_diameter_category)) %>% 
     ggplot(.) +
     geom_bar(aes(x=tree_diameter_category, 
                      fill=tree_health_simplified),
                      stat="count", position = "fill")  +
    labs(x="Categoría diamétrica", y= "Porcentaje de árboles") +
    theme(text = element_text(size = 15)) +
    ggtitle("Estado de salud de los árboles según su diámetro") + 
     theme(plot.title = element_text(lineheight=1.1, face="bold")) + 
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
    facet_grid(. ~ reforested,
    labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
    theme(axis.text.x = element_text(angle=90)) + theme_bw()

Ahora evaluaremos el nivel del daño por ozono, es decir, qué porcentaje del árbol se encuentra dañado, dentro del subset de árboles que están dañados.

## base data
p_od<- muestreo_tidy %>% filter(!is.na(ozone_damage_percentage)) %>%
            ggplot() +
            scale_fill_manual(values= my_cols2, 
                              breaks = desired_order_percentage,
                              labels = c("menos de 10%", "10 a 40%", "40 a 50%",
                                         "50 a 70%", "más de 70%"),
                              name= "Porcentaje del árbol \n dañado por ozono") +
            theme_bw() + theme(text = element_text(size = 20))

Las siguientes figuras muestran que los árboles más viejos tienden a tener mayor porcentaje del árbol con daño por ozono, aunque los árboles más dañados (>70% del árbol) se encuentran en edades jóvenes e intermedias.

p_od +
  geom_bar(aes(x=tree_diameter_category,
               fill=ozone_damage_percentage)) +
  labs(x="Categoría diamétrica", y= "Número de árboles") +
  ggtitle("Porcentaje del árbol dañado en árboles con daño por ozono") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))

Viéndolo por la edad de los árboles (en árboles <15 años):

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage)) +
  labs(x="Número de nodos", y= "Número de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))
## Warning: Removed 100 rows containing non-finite values (stat_count).

En términos porcentuales:

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage),
           position = "fill") +
  labs(x="Número de nodos", y= "Porcentaje de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))
## Warning: Removed 100 rows containing non-finite values (stat_count).

Al explorar la distribución entre árboles de regeneración natural y reforestaciones notamos una tendencia de mayor porcentaje de daño en las plantas reforestadas en comparación a plantas de regeneración natural de la misma edad:

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage),
           position="fill") +
  labs(x="Número de nodos", y= "Porcentaje de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))
## Warning: Removed 100 rows containing non-finite values (stat_count).

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage)) +
  labs(x="Número de nodos", y= "Porcentaje de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural")))
## Warning: Removed 100 rows containing non-finite values (stat_count).

Sin embargo, parece que la exposición (una planta se consideró expuesta si no estaba a la sombra inmediata de otro árbol o arbusto) puede influir en qué tan dañada está una planta:

# plot
p_od +
  geom_bar(aes(x=tree_exposition,
               fill=ozone_damage_percentage),
           position = "fill") +
    facet_grid(. ~ reforested, 
             labeller = as_labeller(c("yes" = "reforestado",
                                    "no" = "natural"))) +
  scale_x_discrete(breaks=c("cover", "exposed"),
        labels=c("cubierto", "expuesto")) +
    labs(x="", y= "Porcentaje de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold"))

Finalmente, hay una tendencia a que mayor porcentaje del árbol se encuentre dañado cuando además del ozono, el árbol presenta síntomas de alguna otra fuente de estrés:

p_od +
  geom_bar(aes(x=tree_nodes,
               fill=ozone_damage_percentage),
           position="fill") +
  labs(x="Número de nodos", y= "Porcentaje de árboles") +
  ggtitle("Porcentaje del árbol dañado por ozono en árboles <15 m") + 
  theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  facet_grid(. ~ tree_health_simplified, 
             labeller = as_labeller(c("ozone" = "Sólo daño por ozono",
                                    "ozone_and_other" = "Daño por ozono y otro")))
## Warning: Removed 100 rows containing non-finite values (stat_count).

Distribución del estado de salud en diferentes reforestaciones.

refos<-filter(muestreo_tidy, reforested=="yes")
nrow(refos)
## [1] 342

En esta seccion exploramos el estado de salud de las plantas concentrándonos únicamente en las plantas reforestadas, que en total son 342.

refos %>% 
     ggplot(.) +
     geom_bar(aes(x=reforestation_year, 
                      fill=tree_health_simplified),
                      stat="count")  +
    scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud del árbol") +
    theme_bw() +
    labs(x="Etiqueta de reforestación", y= "Número de árboles") +
    ggtitle("Estado de salud de los árboles según la reforestación") +
    theme(plot.title = element_text(lineheight=1.1, face="bold")) + 
    theme(text = element_text(size = 25))